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SUMMARY 


We present ideas on how to use wavelets in the solution of boundary value ordinary differential 
equations. Rather than using classical wavelets, we adapt their construction so that they become 
(bi) orthogonal with respect to the inner product defined by the operator. The stiffness matrix in a 
Galerkin method then becomes diagonal and can thus be trivially inverted. We show how one can 
construct an O(N) algorithm for various constant and variable coefficient operators. 

INTRODUCTION 


The purpose of this paper is to use wavelets in the solution of certain linear ordinary differential 
equations of the form 

771 

Lu(x) = f(x ) for x G [0, 1], where L — y^a,(x) P J , 

i= o 


and with appropriate boundary conditions on u(x) for x = 0, 1. 

Currently there exist two major solution techniques. First, if the coefficients a.j(x) of the 
operator are constants, then the Fourier transform is well suited for solving these equations. The 
underlying reason is that the complex exponentials are eigenfunctions of a constant coefficient 
operator and they form an orthogonal system. As a result the operator becomes diagonal in the 
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Fourier basis and can thus trivially be inverted. The numerical algorithm then boils down to 
calculating the discrete Fourier transform of the right hand side, dividing each coefficient by its 
corresponding entry in a diagonal matrix and finally taking the inverse Fourier transform to obtain 
the solution. This can be done quickly using the fast Fourier transform which has a complexity of 
N log N , where N is the number of unknowns in the discretization. 

If the coefficients are not constant one typically uses finite element or finite difference methods to 
discretize the problem. We focus here on finite element methods. Define the opr mini- inner product 
associated with an operator L by 

((u,v)) = { Lu,v ). 

A weak solution u can be found with a Petrov-Galerkin method, i.e. consider two spaces S and <S* 
and look for a solution u G S such that 


((u,v)) = if,v), 

for all v in S*. If S and S’ are finite dimensional spaces with the same dimension, this leads to a 
linear system of equations. The matrix of this system, also referred to as the *li.JJnr$$ matrix, has as 
elements the operator inner products of the basis functions of S and S*. 

Traditionally one uses very local finite elements such that the stiffness matrix has a banded 
structure. The linear system can then be solved efficiently with an iterative method. These classical 
finite elements however have the disadvantage that the stiffness matrix becomes ill conditioned as 
the problem size grows. This slows down the convergence speed of the iterative algorithm 
dramatically. It is well understood by now that this can be solved with multiresolution techniques 
such as multigrid or hierarchical basis functions [1, 2]. Multiresolution finite element bases can 
provide preconditioners which result in a uniformly bounded condition number, see e.g. [3, 4, 5]. 
The convergence of the linear system is then independent of the problem size. 

The research presented here is motived by the question of how good wavelets are for the solution 
of ordinary differential equations. We know that there are basically four main properties of 
wavelets; namely, they provide a multiresolution basis for a wide variety of function spaces, they are 
local in both space and frequency, they satisfy (bi) orthogonality conditions and fast transform 
algorithms are available. Because of these properties, wavelets have already proven to be a valuable 
substitute for the Fourier transform in many applications. 

One possible idea, as proposed by several researchers, is to use wavelets as basis functions in a 
Galerkin method. This has proven to work and results in a linear system that is sparse because of 
the compact support of the wavelets, and that, after preconditioning, has a condition number 
independent of problem size because of the multiresolution structure. However, in this setting the 
wavelets do not provide significantly better results than more general multiresolution techniques 
(cfr. supra) and in fact one of their major properties, namely their (bi)orthogonality, is not 
exploited at all. 

Three questions are addressed in this research. The first, how can one make use of the 
(bi)orthogonality property of the wavelets? The second, which operators can be diagonalized by 
wavelets? The last, are fast algorithms available and what is their complexity? 
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PRELIMINARIES 


Notation and definitions 

Much of the notation will be presented as we go along. Here we just note that the inner product 
of two square integrable functions /, <7 € L 2 (IR ) is defined by 

(f,g) = J ^ f{x)g{x)dx, 

and that the Fourier transform of a function / is defined as 

f(u) = J + " f(x)e~ iux dx. 

We say that a function w is an L-spline if 

L“Lw = 0 and w £ C 2m ~ 2 , 

where L* is the adjoint of L, a linear differential operator of order m. This definition leads to the 
classical piecewise polynomial splines in case L = D m . 

Multiresolution analysis 

We give a brief review of wavelets and multiresolution analysis. For more information one can 
consult [6, 7, 8, 9]. A multiresolution analysis of L 2 (IR) is defined as a set of closed subspaces Vj, 
with j £ that exhibit the following properties: 

1. Vj C V j+1 , 

2. v(x) e V, v(2x) £ Vj+ 1 and v(x) £ Vo v{x + 1) £ Vo, 

+ X' +CXL 

3- U Vjis dense in L 2 (IR ) and Vj = {0}, 

j=-oc, 

4. A scaling fund ion <f)(x) £ V Q exists such that the set of functions {(pjj(x) \ l £ ZZ}, with 
<f>jj(x) = \j2i 0(2- 7 x — 1), is a Riesz basis of Vj. 

As a result there is a sequence {hk \ k £ ZZ) such that the scaling function satisfies a refinement 
equation 

4>{x) = 2Y J h k <t>{2x - k). (1) 

k 
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Define W } now as a complementary space of V 3 in Vj+\, such that V }+ \ = Vj © Wj (® stands for 
direct sum) and, consequently, 

©Wj = L\IR). 

J 

Note that this definition of Wj as a complementary space is non unique. 

A function ip(x) is a -warn let if the set of functions {ip(x — l) \ l E 2Z} is a Riesz basis of W 0 . The 
set of wavelet functions {ifij } i(x) \ l, j E is then a Riesz basis of L 2 (/jR). Since the wavelet is an 
element of V\, it too satisfies a refinement relation, 

ip(x) = 2^2g k (j>{2x - k). (2) 

k 

There are dual functions 4>j,i(x) = y/^4>(2 j x - l ) and ^pj,i{x) = - l) that exist so that the 

projection operators Pj and Q 3 onto Vj and Wj, respectively, are given by 

Pjf(x) = and Qjf( x ) = 

i i 

The basis functions and dual functions are biorthogonal, 

= S i-i' and = Sj-i'&i-v- (3) 

If the basis functions are orthogonal, they coincide with the dual functions and the projections are 
orthogonal. 

The dual scaling function and wavelet satisfy 

4>{x) = 2^2, h k 4>(2x — k), 4>{x) = 2^2 9k <f>(2x — k), (4) 

k k 

and 

<f>(2x - k) = h k- 2 l 4>i x - 0 + 9k-2l ^{x - /). (5) 

l l 

Taking the Fourier transform of the refinement equations (1) and (2) yields 

4>{u) = h(u>/2)4>(u)/2) with h(ui ) = h k e~ ikw 

k 

and 

V»(w) = g(u>/2)ip(uj/2), with g( u>) = Y^ 9 ke~ iku . 

k 

Here h{u) and g(ui) are 27r-periodic functions that correspond to discrete filters. Similar definitions 
and equations hold for the dual functions. A necessary condition for biorthogonality is then 

Vu G IR : m(u;)m f (u;) = 1, 

where 

m(u>) = 

and similarly for rh(u). The existence of the dual filters is guaranteed by the following lemma: 
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Lemma 1 The spare generated by the set of functions I l € com pit nit nls Vj in Vj+i if and 
only if 6(u) = <l>. I m(u) does not vanish. 

The following statements are now equivalent : 

• The dual wavelet has M vanishing moments. 

• Any polynomial with degree less than M can be written as a linear combination of the 
functions <f>j,i(x) with l € ZZ. 

• If f E C M , then the error of the approximation ||/ — Pjf\\ decays as 0(h M ) with h = 2~i . 
These statements are also equivalent with the Strang-Fix condition [10]. 


The fast wavelet transform 


Since Vj is equal to i © W )_ lt a function Vj 6 Vj can be written uniquely as the sum of a 
function Vj_i 6 V)_i and a function Wj_\ € Wj-\'. 

Vj( x ) = Jl^j,k<f>j,k{ x ) = v i-i( x ) +w j - l (x) 

k 

= 53 <f>j- iA x ) + 53 »j-hi i’j^A x )- 

i i 

There is a one-to-one relationship between the coefficients in the different representations. The 
decomposition formulae can be found using (4): 

u j-i,t — V2 hk -21 u j,k , and fij-ij = v5 9 k- 2 i Vj,k- 

k k 

The reconstruction step involves calculating the Vj^ from the Vj-\,i and the Using (5) we have 

u j,k = A2 ^ hk -21 + v5 9k-2i 

i i 

When applied recursively, these formulae define a transformation, the fast wavelet transform [8, 11]. 
The decomposition step consists of applying a low-pass ( h ) and a band-pass (g) filter followed by 
downsampling (i.e. retaining only the even index samples). The reconstruction consists of 
upsampling (i.e. adding a zero between every two samples) followed by filtering and addition. Note 
that the filter coefficients of the fast wavelet transform are given by the coefficients of the 
refinement equations. 

There are many constructions of wavelets. Here we shall only consider compactly supported 
wavelets as in [12, 13]. In this case the filters used in the fast wavelet transform are finite impulse 
response filters and a fast accurate implementation is assured. 
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General idea 


We shall assume that L is self-adjoint and positive definite and, in particular, we can write 

L = V*V, 

where V* is the adjoint of V. We call V the square root operator of L. Suppose that { } and 
for an appropriate range of indices, are bases for S and S* respectively. The entries of the 
stiffness matrix are then given by 




Now, the idea is to let 

= V- 1 ^ and = V 1 ^ , 

where ip and ip are the wavelets of a classical multi resolution analysis. Because of the 
biorthogonality (3), the stiffness matrix becomes a diagonal matrix which can trivially be inverted. 
This avoids the use of an iterative algorithm. We will call the \1> and 'f r * functions the operator 
wavelets and the ip functions the original wavelets. The operator wavelets are biorthogonal with 
respect to the operator inner product, a property we refer to as operator biorthogonal. 

This idea can be powerful, but there are a few problems. First of all one has to check whether 
the operator wavelets still provide an multiresolution analysis where the successive approximations 
to a general function converge sufficiently fast (efr the Strang-Fix condition). Secondly one has to 
construct a fast wavelet transform for this operator multiresolution analysis. We want operator 
wavelets to be compactly supported and to be able to construct compactly supported operator 
scaling functions We will see that the latter is not as simple as just applying V~ l to the 
original scaling functions. 

The analysis is relatively straightforward for simple constant coefficient operators such as the 
Laplace and polyharmonic operator. For more general constant coefficient operators, we will show 
that one needs to modify the construction of the original wavelets for the operator wavelets to 
satisfy all the desired properties. We will discuss the Helmholz operator as a typical example. At 
the end of the paper we shall consider a variable coefficient operator. 

A similar idea was described in [14, 15]. However there only the operator wavelets of different 
levels are operator orthogonal and not the ones from the same level. As a result, one does not 
obtain a full diagonalization, but rather a decoupling of equations corresponding to different levels. 

Our idea is different from the technique presented in [16]. There wavelets are used to efficiently 
compute the inverse of the matrix that comes from a finite difference discretization. It is also shown 
that the wavelets provide a diagonal preconditioner which yields uniformly bounded condition 
numbers. 

In [17, 18] antiderivates of wavelets are used in a Galerkin method. This parallels our 
construction in the case of the Laplace or polyharmonic operator. 
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LAPLACE OPERATOR 


The one dimensional Laplace operator and its square root axe 

L = -D 2 and V = D. 

The associated operator inner product is therefore ({ u,v )) = (u',v') . Since the action of V is 
simply taking the antiderivative here, we define the operator wavelets as 

V(x) = [ X ^(t)dt, and ®*(x) = [ rp{t)dt. 

7—00 7-co 

The operator wavelets are compactly supported because the integral of the original wavelets has to 
vanish. Also translation and dilation invariance is preserved, so we define 

= ^(2 j x-l) and ^- >{ (a;) = \t*(2 j x-l). 

It is then easy to see that 

{{ )) = 2 j 6j_j> 6i-y for j,j\ 1,1' e2Z. 

This means that the stiffness matrix is diagonal with powers of 2 on its diagonal. 

We now need to find an operator scaling function $. The antiderivative of the original scaling 
function is not compactly supported and hence not suited. We instead construct the operator 
scaling function 4> by taking the convolution of the original scaling function with the indicator 
function on [0, 1], 

$ = <j>*X[ o.ip 

and similarly for the dual functions. We will show that these functions indeed generate a 
multiresolution analysis. To this end define 

Vj = clos span {$j t i | l € S} and Wj = clos span {^j,i \ l G 2Z}. 

We show that the V 3 spaces are nested and that Wj complements Vj in Vj+\. 

In the Fourier domain we have 

1 — p~ iw ~ ~ 1 - 

$(oj) = Moj) and ^(u;) = 

iW VJ 

A simple calculation shows that the operator scaling function satisfies a refinement equation 

l+e _iw 

$(o,) = $(u/2) H(u/2) with H(u > ) = h(u>). 

Consequently, the Vj spaces are nested. If we can find a function G such that 

*(u) = $(u/2)G(u)/2), 
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then this implies that Wj is a subset of V j+1 . It is easy to see that this holds with 


G(u) 


1 

2{1 - e~ iw ) 




This function is well defined because g(0) = 0. 
The space Wj complements V 3 in Vj+i if 


A (uj) = det 


H(u) H(u + tt) 
G(u) G( u/ + 7r) 


does not vanish, see lemma 1. In fact, we readily see that A(u;) = 6(u)/ 4, and this cannot vanish 
since (f> and ip generate a multiresolution analysis. The construction of the dual functions and 
from 4> and ip is competely similar. The coefficients of the trigonometric functions H, H' , G and 
G* now define a fast wavelet transform. 

Note that there is no reason why the operator scaling functions should be operator biorthogonal 
and in fact one can prove that this never happens. Note also that if true, this property would make 
the use of wavelets superfluous. 


Algorithm 

We will describe the algorithm in the case of periodic boundary conditions. This implies that the 
basis functions on the interval [0, 1] are just the periodization of the basis functions on the real line. 
Let S = V n and consider the basis {<t> ni | 0 ^ l < 2 n }. Define vectors b and x such that 

b i = ( /. K,i ) » and u = Y, x i 

(=0 

The Galerkin method with this basis then yields a system 

Ax = b with A k>l = <($ ni ,,$ ni *)). 

As we mentioned earlier, the matrix A cannot be diagonal. Also its condition number grows as 
G( 2 2n ). Consider now the decomposition 

v n = Vo© wbe-’-ew'fi-i, 

and the corresponding wavelet basis. The space V 0 has dimension one and contains constant 
functions. We now switch to a one index notation such that the sets 

{1, V jtl \0^j<n,0^l<2 j } and {** | 0 ^ A; < 2 n } 

coincide. Define the vectors b and x such that 

2 " — 1 

h - (/, ^/) and u = x t 

t=o 

We know that there exists matrices T and T* such that 

b = T* b and x — Tx 


266 


The matrix T' corresponds to the fast wavelet transform decomposition with filters H' and G" and 
T corresponds to reconstruction with filters H and G. The complexity of the matrix vector 
multiplication is O(N), N = 2 n . In the wavelet basis the system becomes 

Ax = b with A = T* AT and Akj = ^n,k )) • 

Since A is diagonal, it can be trivially inverted and the solution is then given by 

a; = TA~ x T'b. 

This means that one has to calculate the wavelet decomposition of the right hand side, divide each 
coefficient by its corresponding diagonal element and reconstruct to find the solution. The 
complexity is O(N). 

The constant basis function of Vo has a zero as diagonal element and its coefficient is thus 
undetermined. Note that this leads to an inconsistency if the integral of / does not vanish. 


Boundary conditions 


Our general idea to deal with boundary conditions is to let the operator wavelets satisfy the 
homogeneous boundary conditions and to let the component in the Vo space satisfy the imposed 
boundary conditions. This requires the use of special boundary wavelets as described in [19]. With 
only a slight change of basis one can then incorporate Dirichlet, Neumann, mixed and periodic 
boundary conditions. The details of this construction go beyond the scope of this paper. We will 
describe the construction in some specific cases. 


Example 


In this section we shall take a look at a simple example, namely the basis we get by starting 
from the Haar multiresolution analysis, where 

<f> = X[o,i] and ^{x) = <f>(2x) - <p{2x - 1). 


Define the hut function as 

A = X[o,i] * X[o,i] ■» such that $ = A and \V(x) = A(2x). 

The original wavelets are orthogonal and as a consequence the basis functions and dual functions 
coincide. 

The operator scaling functions can represent linears which means they satisfy the Stang-Fix 
condition with M = 2 and the convergence is of order h 2 . One can prove that higher order wavelets 
with more vanishing moments (M) will in general not yield faster convergence because the solution 
u is not smooth enough. The underlying reason is that the solution u belongs to the Sobolev space 
W 2 . One can get faster convergence only by imposing extra regularity conditions on the right hand 
side. So in a way this basis seems to be the most natural one to work with. Note that these 
piecewise linear basis functions are local solutions of the homogeneous equation such that the 
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Figure 1: Basis for Dirichlet problem. 


Figure 2: Basis for Neumann problem. 


operator seeding functions and wavelets are V-splines. This basis also coincides with Yserentant’s 
hierarchical basis. 

Figure 1 shows the basis functions in the case of Dirichlet boundary conditions and n = 3. The 
left part are the bases for the spaces V 0 up to V 3 while the right part are the bases for W 0 up to W 2 , 
which provide the diagonadization. The coefficients of the two functions in the Vo space are 
determined by the boundary conditions. The fast wavelet transform differs from the periodic 
algorithm here in the sense that different coefficients are used for the wavelets at the boundary. 
Note the “hedf hat” functions here. The basis in case of the Neumann problem is shown in figure 2. 
The boundary conditions are handled by the two functions in the Vi space. Again the coefficient of 
the constant is undetermined. The algorithm leads to an inconsistency in case the integral of / is 
not equal to it'(l) — u'(0). Note that in both cases the operator wavelets satisfy the homogeneous 
boundary conditions. 

MORE GENERAL CONSTANT COEFFICIENT OPERATORS 


The polyharmonic operator 


The polyharmonic equation is defined as 


-u (2m) = /, 

and the square root operator is now V = D m . The operator scaling function $ is now m times the 
convolution of the original scaling function (p with the box function and the operator wavelet 'k is 
m times the antiderivative of the original wavelet ip. In order to get a compactly supported wavelet, 
the original wavelet now needs to have at least m vanishing moments, a property which can be 
satisfied by all known wavelet families. The construction and algorithm are then completely similar 



Figure 3: The refinement relation for the piecewise exponentials. 


to the case of the Laplace operator. 


The Helmholz operator 

The general definition of the one dimensional Helmholz operator is: 

L = —D 2 + k 2 such that V = D + k. 

Here we shall assume that k = 1 which can always be obtained from a simple transformation. 
Observe that V = D + I = e x De x and thus V~ l = e~ x D~ l e x . One easily verifies that applying 
V- 1 to a wavelet will not necessarily yield a compactly supported function since e x i)j t i in general 
does not have a vanishing integral. Therefore we let 'Fj,; = V~ l e~ x ipj t i = e ~ x If ipjj has a 
vanishing integral, then is compactly supported. 

In order to diagonalize the stiffness matr ix, the original wavelets now need to be orthogonal with 
respect to a weighted inner product with weight function e~ 2x because 

(( j, i, )} = J + " e~ 2x dx. 

Finding such wavelets is a hard problem to solve in general. Inspired by the Haar basis, we 
construct a solution where the orthogonality of the wavelets on each level immediately follows from 
their disjoint support, by letting supp^^ = [2" J Z, 2~*(l + 1)]. To get orthogonality between the 
different levels, we need that Vj is orthogonal to Wj> for f ^ j or 

J e~ 2x <pj,t(x) ipj' t i'(x) dx = 0 for j’ ^ j. 

We now let the scaling function coincide with e 2x on the support of the finer scale wavelets, 

j 2x 

= e Xj.i , 
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where Xj,i is the indicator function on the interval [2 j l, 2 J (l + 1)], normalized such that the 
integral of the scaling functions is a constant. As in the Haar case we choose the wavelets as 


V’j.i = <t>j+l,2l - <j>j+l,2l+l, 

so that they have a vanishing integral. The orthogonality between levels now follows from the fact 
that the scaling functions coincide with e 2x on the support of the finer scale wavelets, and from the 
vanishing integral of the wavelets 

/ -fOO _ f+O O _ /*+ OC* _ 

e_2;r 0j,i( x ) fpj>,i'(x)dx = / = / ip r<l '(x) dx = 0. 

-OO J-o O J—oo 

One can see that the operator wavelets are now piecewise hyperbolic functions (piecewise 
combinations of e x and e~ T ). The scaling functions are chosen as 

D~ l — (f>jj + 1 ) so that 


With the right normalization, one gets 


*jA x ) 


sinh(x — 12 J ) 
sinh(2~ J ) 

i sinh((f + 2)2 _J — x) 
sinh(2 - '?) 


for x € [I2~ j , ( l + l)2~ j ] 
for i£[(l + l)2- j , (l + 2)2~ j \ 


0 


elsewhere. 


The operator scaling functions on one level are translates of each other but the ones on different 
levels are no longer dilates of each other. They are supported on exactly the same sets as the ones in 
figure 1 and they roughly look similar. The operator scaling functions satisfy a refinement relation 

*u = 

k—0 


with 

Hi = = sinh(2 _ - 7_1 )/sinh(2 _ - 7 ) and H{ = 1. 

Figure 3 shows the refinement relation for the scaling functions. The 3 finer scale functions are not 
the dilates of the coarse scale one but they still add up to it. 

The Helmholz operator in this basis of hyperbolic wavelets again is diagonal and the algorithm is 
completely similar to the Laplace case. The only difference in implementation is that the filter 
coefficients H 3 k used in the fast wavelet transform now depend on the level. 

Note that these functions again are F-splines and, in a way, are the most natural to work with. 
Also note that 

lim 2 -J x) = A(x). 

j—x 

Despite the fact that the Strang-Fix conditions are not satisfied, one can prove that the 
convergence is still of order h 2 . 

So we can conclude that a wavelet transform can diagonalize constant coefficient operators 
similar to the Fourier transform. The resulting algorithm is a little faster (O(N) instead of 
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C?(A T log N)). This gain in speed is a consequence of the subsampling of the coarser levels in the 
wavelet transform (the ones that correspond to the low frequency components of the solution) 
which is not present in the Fourier transform. Also boundary conditions are taken care of more 
easily than in the Fourier case. 


VARIABLE COEFFICIENTS 


Naturally, the next question is how to use wavelets for variable coefficient operators. The 
underlying reason why wavelets can diagonalize constant coefficient operators is their locality in the 
frequency domain. We want to understand if we can exploit the localization in space to diagonalize 
variable coefficient operators. The answer is (perhaps quite surprisingly) yes and this really justifies 
the use of wavelets for differential equations. No other technique (to our knowledge) has been able 
to accomplish this. 

We take a closer look at the following operator 

L = —D p 2 (x)D, 

where p is sufficiently smooth and positive. The square root is now V = pD and V -1 = D~ l 1/p. 
The rest of the analysis is very similar to the case of the Helmholz operator. Applying V -1 directly 
to a wavelet does not yield a compactly supported function. We therefore take pipjj 

which implies that the wavelets need to be (bi) orthogonal with respect to a weighted inner product 
with p 2 as weight function. We use the same trick as for the Helmholz equation to construct such 
functions. This means that we let the scaling functions <f> 3 j coincide with 1/p 2 on the dyadic 
interval [2~ j l,2~ j (l + 1)] and normalize them such that they have a constant integral. We then take 
the wavelets ip jit to be equal to ~ <f>j+ i,2t+i so they have a vanishing integral and the operator 

wavelets are compactly supported. The operator wavelets are now piecewise functions that locally 
look like AP + B where P is the antiderivative of 1/p 2 and again are V-splines. Their support also 
coincides with the support of the functions of figure 1, and since p is smooth they will converge to 
hat functions as the level goes to infinity. The operator wavelets are neither dilates nor translates of 
one function, since their behavior locally depends on p. This is not a problem because they still 
generate a multiresolution analysis and satisfy refinement relations. The coefficients in the fast 
wavelet transform are now different everywhere and they depend in a very simple way on the Haar 
wavelet transform of 1/p 2 . The entries of the diagonal stiffness matrix can be calculated from the 
wavelet transform of 1/p 2 - The algorithm is completely similar to previous cases and is of order N. 
Boundary conditions are as easy to handle as in the case of the Laplace operator. Note that the 
operator scaling functions do not satisfy the Strang-Fix conditions. It is however again possible to 
prove that the method has a convergence of order h 2 . As mentioned earlier, higher convergence 
orders can not be obtained in general. 

NUMERICAL EXAMPLE 


We solve the equation 

—D ef Du(x ) = e x> (sin(x)(3x 2 - 2) + cos(x)(2x - 2x 3 )) / x 3 , with u(0) = 1 and u(l) = sin(l), 
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l 

L^, error 

1 

1.22e-02 

2 

3.37e-03 

3 

8.66e-04 

4 

2.18e-04 

5 

5.45e-05 

6 

1.36e-05 

7 

3.41e-06 

8 

8.52e-07 

9 

2.13e-07 


such that the exact solution is given by u(x) = sin(x)/x. The error of the numerically 
computed solution is a function of the number of levels (l) shown in the above table . Each time 
the number of levels is increased the error is divided almost exactly by a factor of 4, which agrees 
with the G(h 2 ) convergence. 


CONCLUSION 

In this paper we showed how wavelets can be adapted to be useful in the solution of differential 
equations. Like the Fourier transform, wavelets can diagonalize constant coefficient operators. The 
resulting algorithm is slightly faster. The main result however is that even non-constant coefficient 
operators can be diagonalized with the right choice of basis which evidently yields a much faster 
algorithm than more classical iterative methods. 

This technique can also be applied to the solution of implicit time stepping discretizations of 
equations of the form du/dt = Lu -t- f even when L is non-linear. Future research includes the 
study of non self adjoint operators where a splitting L = VV“ is needed and the study of the 
possible generalization of these ideas to partial differential equations. 
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